c^ 



Abstract 



"(Si". A unified framework for the computation of 

O : polynomial quadrature weights and errors* 

c^ I Mario M. Graca and M. Esmeralda Sousa-Dias ^ 

March 1, 2013 
O 

m 

< 

^^ ' For the class of polynomial quadrature rules we show that conveniently 

rS^ [ chosen bases allow to compute both the weights and the theoretical 

j^ ' error expression of a n-point rule via the undetermined coefhcients 

method. As an illustration, the framework is applied to some classical 
rules such as Newton-Cotes, Adams-Bashforth, Adams-Moulton and 
Gaussian rules. 

(N 
> 

lO ' Key-words: Quadrature, undetermined coefficients method, degree of 

^^ . precision, orthogonal polynomial. 

'Nh ■ MSC2010: 65D30, 65D32, 65D05, 42C05. 

cn 

o 

CN ■ 1 Introduction 

Most numerical integration schemes widely used in the applications are 
►v> ' based on the class of polynomial quadrature. Here we present an unified 

^ . framework for the simultaneous computation of weights and theoretical er- 

ror of a n-point polynomial quadrature rule Qnif)- We aim to bring together 
the computation of the pair (weights, error) = {Wn, En) under a simple and 
unified framework. Our approach shows that it is possible to choose a par- 
ticular polynomial basis relative to which the weights Wn are the solution 
of an upper triangular linear system, enabling so their recursive computa- 
tion. At the same time, by extending conveniently the referred basis, we 
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obtain a (computable) criterium for the degree of precision of the rule and 
consequently the respective error expression En. 

The upper triangular system, whose solution is the vector Wn of weights, 
is obtained by the technique known in numerical analysis as undetermined 
coefficient method (not to be confused with the same named method for 
ODE's). The undetermined coefficients method (UCM) is often employed 
for weights's computation, however this method has been set aside for the 
determination of En (see for instance the remarks in Gautschi [7], p. 176, 
and Evans |5], p. 69). To the best of our knowledge, our scheme for the 
simultaneous computation of the weights and error via the UCM is new. 

In numerical quadrature one aims to approximate an integral /(/) = 
J^ f{x)dx by a n-point rule Qn{f) = X]fc=i '^^i/l^i); ^^d to obtain an ex- 
pression or a bound for the error En{f) = I{f) — Qnif)- The coefficients Oj 
in Qnif) are called the (quadrature) weights and the points Xj in a given set 
{xi,...,x„} are usually called nodes (a.k.a. quadrature points, or abscis- 
sae). The rule Qnif) is said to be polynomial interpolatory when it is exact 
for polynomials of degree less or equal to n — 1. That is, when Qnip) = Hp) 
for all polynomials p of degree less or equal to n — 1. Hereafter when we 
refer to a rule we mean a polynomial quadrature rule. 

Quadrature has a long history and is a classical subject, although it 
is still nowadays an active area of research (general references are, for in- 
stance, Atkinson [1], Davis and Rabinowitz [3], Gautschi [7] and Krylov |10)). 
There are a priori two main goals when dealing with numerical quadrature: 
the determination of the weights and respective theoretical error expression. 
Traditionally, in order to obtain the weights, one first determines the interpo- 
lating polynomial for a given panel {(xi, /(xi)), (x2, fix2)) • • • , (x„, /(x„))} 
and afterwards this polynomial is integrated. The interpolating polynomial 
is unique (and so are the weights of the rule) but it may be written in several 
different forms named after Lagrange, Newton, Hermite, Chebyshev, etc. Of 
course, this means that we can always write the interpolating polynomial in 
distinct polynomial bases. 

In the classical literature the theoretical quadrature error is usually de- 
duced from a particular expression representing the interpolation error - 
see Steffenson [14] for a survey on interpolation theory, and Berezin and 
Zhidkov [2], Ch. 2, for its applications. Consequently, the quadrature error 
deduction is based upon a panoply of analytical tools such as mean value 
theorems, integration by parts, Peano kernels, Euler-MacLaurin formula and 
so on. The usual approach for the computation of the pair (W„, £'„), where 
Wn = (oi, . . . , a„), is done on a rule by rule basis leading usually to lengthy 



and cumbersome calculations even for a small number n of nodes. 

At the heart of our approach are conveniently chosen bases of polynomi- 
als relative to which we apply the undetermined coefficients method in order 
to compute the pair {Wn,En)- Only elementar results from linear algebra 
and numerical analysis are needed. 

The paper is organized as follows. There are two sections, in the first 
one we prove the main result (Theorem 12. ip . which points out an algorithm 
for computing the pair (Wn,En). In the following section we illustrate its 
applicability to well known quadrature rules named after their creators: 
Newton-Cotes, Adams-Basforth, Adams-Moulton, and Gaussian rules. This 
shows the wide applicability of the algorithm previously referred. In par- 
ticular, when applying Theorem 12. II to Gaussian quadratures we recover (in 
Proposition 13. ip a well known closed form for their error expression. 

We barely address here numerical computational issues since they are 
out of the scope of this work, and will be treated elsewhere. However, 
several symbolic/numerical tests were carried out in order to compute pairs 
{Wn,En) for the illustrative quadrature rules mentioned above. For, we 
developed a Mathematica code translating the algorithm in Theorem 12.11 
and we have tested it on all the rules described in Section [3l The respective 
computed weights and error expressions were compared with those in tables 
spread in the literature. In particular, for a number of nodes 2 < n < 256, 
the algorithm has been used to compute 100-digits precision pairs {Wn,En) 
for the ubiquitous and indispensable Gauss-Legendre rule. 

2 Weights and error for a polynomial rule 

In order to approximate the integral /(/) = / f{x)dx, where / is assumed 
to be an integrable function in (a, b), and of class C^{[a, b]) (where the integer 
k will be taken appropriately to each case), one constructs the rule 

Qnif) = aif{xi) + a2f{x2) H h anf{xn), (1) 

from a given set {xi, X2, . . . , x^}, with n > 1. 

The rule Qn{f) is said to have degree of precision d, if I{p) = Qn{p) for 
all polynomials p of degree < d, and I{p) ^ Qn{p) for some polynomial p of 
degree d+1. If d = deg{Qn{f)) denotes the degree of precision of the rule 
Qnif), and / is of class C"^"'"^ in [a,b], the quadrature error is given by 

Enif) = I{f) - Qnif) = ^-^^^^^^^^^^f^'^'\il (2) 



where ^ and Pd+i are respectively a certain point in (a, b) and a polynomial 
of degree d+1 (cf. Gautschi [;7], p. 178). We remark that the degree of 
precision of a rule cannot be greater than (2n — 1) since, if p is a polynomial 
of degree n then /(p^) is positive. 

In numerical analysis the undetermined coefficients method consists in 
solving the linear system arising from the equalities Qnidi) = I{gi), for 
< z < (n — 1), where gi are the elements of a prescribed set of functions. 
In the context of polynomial quadrature this set will be a set of polynomials. 
The next theorem shows that when the UCM is applied to a conveniently 
chosen basis of polynomials, the weights are obtained recursively and, at the 
same time, the degree of precision d of the corresponding quadrature rule is 
easily found. 

We denote by P^ the linear space of real polynomials of degree less or 
equal to s. Suppose that xi,X2, . . . ,x„ are n > 1 distinct real points and 
consider the basis Bq = {ipQ{x), (pi{x), . . . , c^„_i(x)} of Pn-i, given by 

(Pj{x) = ipj-i{x){x — Xj), 1 < J < n — 1. 

The basis Bq is known in some literature as Newton's basis (see Stef- 
fensen [14], p. 23). We complete the basis Bq with other polynomials 
qn{x),qn+i{x), . . . , qn+k{x) to form a basis P„+fc. Notably, 

qn{x) = C^„_i(x)(x - Xn) ^^. 

qj{x) = qj-i{x){x — Xr), j > n + I with r = j (mod n), 

where A; is a nonnegative integer which will be later specified in accordance 
with the particular quadrature rule under consideration. Note that the 
points Xi (1 < i < n — 1) are roots for all the polynomials qj defined in @. 

Theorem 2.1. Let xi,...,x„ be n > 1 distinct real quadrature points, 
Bq = {ipQ{x), . . . ,ipn-i{x)} the basis o/ P„_i defined in ^, and B = 
Bo U {qnix), ... , qn+kix)} the basis ofFn+k, defined in (JH). 

(a) The undetermined coefficient method applied to the basis Bq deter- 
mines (uniquely) the weights ai of the rule ([T]) for approximating the 
integral I{f) = J f{x)dx. These weights are 

_ I{ipn-l) 

I{'Pi-i)-Y2=i+ifi-i^^k)ak . 1 o 1 

Oj = — , « = n- l,n - 2,... ,1. 

'^i-l{Xi) 



(b) The undetermined coefficient method applied to the basis B, deter- 
m,ines the degree of precision d = deg{Qn{f)) as being the integer 
d > n — 1 for which 

I{qd+i) ^ and lilj) = 0; /'''" ^^^ n < j < d. (6) 

(c) If f is of class C^^^dajb]), the error expression is 

EQM)=Cnf'''^'HO, ^^^/^^" = ^^' (7) 

for some ^ G {o-,h). 

Proof, (a) As the rule Qn{f) is polynomial, applying the UCM to the basis 
Bq of the polynomials ipi of degree i < (n — 1), we have 

Qn{Vi) = H^i), for aU < i < (n - 1). (8) 

The conditions ([8]) are equivalent to the linear system Ax = 6 in the un- 
knowns oi, . . . ,a„, where b = {I {ipo) , I {(pi) , . . . , I {(pn-i))'^ ■ The matrix A 
is a n X n upper triangular matrix whose diagonal entries are 1, (pi{x2), • • •, 
ipn-i{xn)- Furthermore, A is obviously nonsingular since Xj+i is not a root 
of ipi, for 1 < i < n — 1, and so det(^) 7^ 0. The (unique) solution of Ax = b 
is obtained by backward substitution and is given by ([5]). 

(b) Applying the UCM to the basis B, we get an overdetermined linear 
system Ax = b, where x is the vector whose components are the weights, 
and the (n + k) x n matrix A and the vector b are, respectively, of the form 
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b= {I{ipo),...,I{ipn-i),I{qn),---,I{qn+k)) , (9) 

{n+k)xn 



where O denotes the zero matrix. The conditions ([6]) in the statement are 
equivalent to say that d = deg(Qn(/)) is the greatest integer ra + fc for which 
the system Ax = 6 is compatible. 

Let us now show that the degree d of the rule Qn{f) is given by the 
conditions ([6]). Consider p to be a polynomial of degree n -\- k for some 
nonnegative integer k and write p as (unique) linear combination of the 
elements of the basis B of P^+fc. That is, 

n—l n+k 

j=0 j=n 



Then, 

n— 1 n+k n—1 

Qn{p) = X] OijQnifj) + ^ OijQniqj) = ^ ajQni^fj), (10) 

i=0 j=n j=0 

where the last equahty follows from the fact that aU points Xi are roots of 
the polynomials qj. Furthermore, by linearity of the operator /, we have 

n— 1 n+k 

j=0 j=n 

As by dD), Qni^j) = I{Vj)i '^6 conclude from pOj) and ([TT]) that 

n+fc 
/(p)-Qn(p) = E"j^(9^)- (12) 

Therefore the degree of precision d cannot be less than (n — 1). As there 
exists a polynomial ps of degree s < 2n such that I{ps) 7^ 0, this implies 
that I{qs) 7^ for an integer s such that n < s < 2n. It also follows easily 
from (fT^ that Qn{p) = I{p)i for all polynomials p of degree < d, if and only 
if I{qj) = for all n < j < d. So, the conditions ([6]) hold. 

(c) As all the quadrature points are roots of the polynomial qd+i, we 
have Qniqd+i) = 0. Then, it follows from ([2]) that the error's expression is 
given by d?]). 

D 

Remarks 2.1. 1) It becomes clear from the proof of Theorem 12.11 the 
reason why the UCM has been set aside for the computation of the 
degree of a rule. Indeed, if one had applied the UCM to another basis 
B, for instance the canonical basis of polynomials, the matrix A in ([9]) 
would not have the zero block (or equivalently the second term of the 
sum in (jlOp does not vanishes), and so it will be impossible to obtain 
the criterium ([6]) for the degree. 

2) Implicit in the proof of the above theorem is the following algorithm. 
Given the nodes xi, . . . ,Xn and the bases Bq and B; (i) Compute the 
weights Wn using the recursion relation ([5]); (ii) Find the first indice 
j > n such that I{qj) 7^ 0; (iii) d = j — 1 is the degree of the rule; (iv) 
Compute the coefficient c,i in d?]); (v) Output Wn and £'„. 



3) As the degree of precision of a rule cannot be greater that 2n, it is only 
necessary to extend the basis Bq to a basis B of P2n- It also follows 
from the proof of Theorem 12.11 that the polynomials qj only need to 
satisfy the requirement of having all quadrature points as roots. 

4) From Theorem 12. !( we can conclude that the quadrature weights aj of 
a polynomial rule do not depend wether the quadrature points belong 
or not to the integration interval [a, b]. However, the value of the inte- 
grals I{qj) do depend on the points's location relative to the interval 
of integration. In the next section we show how the position of the 
quadrature points may influence the degree of the precision of a rule, 
for instance in the cases of the Newton-Cotes and Adams-Bashforth- 
Moulton rules. 

5) The upper triangular system leading to the computation of the weight's 
vector Wn, has determinant equal to det{A) = YYjZi "^ji^j+i)- ^or n 
large and certain sets of nodes we can have det{A) ~ 0. In this case 
we are facing an ill-conditioned problem which raises challenging com- 
putational issues. 

6) It would be interesting to apply Theorem 12. II to polynomial rules other 
than those treated in Section [3l such as the Clenshaw-Curtis, Patter- 
son, Gauss-Konrad, etc. (see Evans [5j for a survey of other available 
rules, and references therein). 

3 Applications to some classical integration rules 

In this section we apply Theorem 12.11 to some classical rules in order to 
obtain their weights and the respective error's expression. To illustrate 
the wide applicability of our result we begin with rules with equally spaced 
nodes (Newton-Cotes and Adams-Bashforth-Moulton) followed by rules with 
unequally spaced points (Gaussian). 

3.1 Weights and error for Newton-Cotes and Adams-Bash- 
forth-Moulton rules 

The abscissae for both Newton-Cotes (open or closed) rules and Adams- 
Bashforth-Moulton are equally spaced. In the rules of the first type all the 
quadrature points belong to the integration interval [a, b] whilst for the sec- 
ond type there are points outside [a,b]. These particularities of the nodes 
have a decisive influence on the rule's degree. The next Lemma shows how 



the abscissae's distribution relative to the integration interval will affect 
the values of the integrals I{qj) in Theorem 12.11 Although the results in 
Lemma 13.11 are scattered in the literature under somehow different formula- 
tions, we include a proof for the sake of completeness. 

Lemma 3.1. Let yi,y2, ■ ■ ■ ,yn be real points such that yi = yi-i + h, for 
some positive constant h, and qn{x) = {x — yi){x — j/2) ■ ■ ■ {x — yn)- 

i) If yi = a and yn = b, then 

^ 1 if n is odd 

qn{x)dx = < , 

I 7^ U if n IS even. 

a) Ifqn+iix) = qn{x){x-y), withy = yi ory ^ [yi,yn], then j^^^ qn+i{x)dx / 
0. 

Hi) If yn-i = a and yn = b, then J^ qn{x)dx ^ 0. 

Proof. The change of variables 7 : x i-^- t, defined hy x = yi + h [t + ^^), 
mapps 2/1 H- -^, 2/2 1-^ -^ + l,...,yn-i 1-^ ^ - l,yn i-> ^, and 
[a,6]^[-(n-l)/2,(n-l)/2]. 

If n is odd, the integral /^ qn{x)dx is equal to the integral of an odd 
function and so it is zero. That is, 

' qn{x)dx = /i"+l /y^ t{t^ - l)(t2 - 22) . . . ft' - ^^^^] dt = 0. 

If n is even, then J qn{x)dx is equal to the integral of an even function: 

qn{x)dx = /i"+i / \{t^- 1/4) (t^ - 9/4) ■■■(t^- ^"""^^ \ dt / 0. 

For ii), we can write y as y = yi + h (^^^) + r, where r is some nonzero 
constant. Thus, {x — y) = {x — (yi + h (^^)) + r, and so 



Un rUn 

qn+iix)dx = / qnix) 
yi -'yi 



Ti — 1 

X- [yi + h' 



dx — r 

's/i 



ryn 

/ qn{x)da 



Applying (i) and the same change of variables as before, it is easy to con- 
clude: (a) When n is odd, then K2 = and Ki is the integral of an 



even function in [— (n — l)/2, (n — l)/2], and so f^" qn+i{x)dx / 0; (b) 
when n is even, then K2 ^ and Ki is the integral of an odd function in 
[-(n - l)/2, (n - l)/2], thus /^^; qn+i{x)dx ^ 0. 

(iii): When x belongs to the open interval {yn-i,yn) the product (x — 
yi){x — ^2) • • • (x — yn-i) is positive and (x — yn) does not vanishes. There- 
fore, the integrand function qn{x) does not change sign in [y„_i,y„], hence 

Closed Newton-Cotes rules 

In order to approximate /(/) = f^ f{x)dx, the interval [a,b] is divided into 
(n — 1) parts of equal length h = ^^Ef ) S'lid the nodes are: Xi = a + ih, for 
i = 0, 1, . . . , (n — l)/i. By a change of variables, one can consider the n nodes 
to be ti = 0, ^2 = 1; • • • , in = "- — 1- The integral /(/) is 

/(/) = / f{x)dx = h r g{t)dt, 

Ja Jo 

with g{t) = f{a + th). The rule Qn{f) is related to the rule 
Qnig) = big{0) + 625(1) + • • • + bngin - 1) 

by Qnif) = hQn{g)- The elements of the basis B^ in ([3]) are (po{t) = 1, 
(/?j(t) = Lpj^i{t){t — {j — 1)), for J = 1, 2, . . . , n — 1. The polynomials g„ and 
Qn+i in dU are 

g„(t) = t(t - 1) • • • (t - (n - 1)), g„+i(t) = t2(t _ 1) . . . (t _ (n - 1)). 

By ([U of Theorem I2.H the weights 6j for Qn (g) are given recursively by 

liifn-l) _ fo~ (Pn-l{t)dt 



bn 



ipn-i{n- 1) 



(n-1)! 



{^- 1)! 

where we have applied the fact that ipk{k) = k\. The weights for the rule 
Qnif) are a^ = hbi. The degree of Qnig) is at least (n — 1), and from 
Lemma [3T] -(i) the integral I(g„(t)) = Jq i(*~l) ' ' ' it — in — l))dt vanishes 
when is odd and is nonzero when n is even. Also, by (ii) of the same Lemma, 



I{qn+i{t)) 7^ 0. Then, by Theorem 1 2. H the degree d and the error expression 
En{f) for a n-point Newton-Cotes rule are respectively, 

{d = n — 1, if n is even 
d = n, if n is odd 

and 

EM) = Hf) - QnU) = h{I{g) - Qn{g)) 

~^{d + i)i^ ^^'~ {d + iy/ ^^^' 

where the last equality follows from the fact that g^^'{t) = h^f^{a + th), 
and we are assuming that / G C'^~^^{[a,b]). That is, 

EM) = h^^'cr,f(^^'HO, withc„, = ^^. (13) 

Using an analogous procedure, it is now a simple exercise to obtain the pair 
(Wn,En) for a n-point open Newton-Cotes rule. 

The Adams-Bashforth-Moulton rules 

The Adams-Bashforth (AB) rules are the core of explicit methods for ODEs 
and the Adams- Moulton (AM) rules play a central role as implicit methods 
for ODEs with the same name (Henrici [9], Gautschi [7], Cheney [3]). 

The (equally spaced) nodes for the AB rule are less or equal to a (where 
[a, 6] is the integration interval), while in the AM case the points a and b are 
quadrature points and there are no other abscissae in the interior of [a, b] . 
The Lemma 13.11 and Theorem 12.11 imply that both the AB and AM rules 
have degree of precision d = (n — 1). 

We now describe in detail the polynomial bases to be used and we de- 
termine the theoretical error expression for both rules. 

The Adams-Bashforth rule 

Given a constant step h and a fixed real number r, consider n > 2 nodes 
xi = r, X2 = T — h, . . . , Xn = {t — {n — l)h). The Adams-Bashforth rule 
aims to approximate the integral /(/) = J^ f{x)dx by 

Qnif) = aif{xi) -\ h a„f{xn)- 
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Making the change of variables defined hy x = r -\- ht^ we have 

1(f) = / f{x)dx = h f{T + ht)dt = h I g{t)dt. 
Jt Jo Jo 

The nodes Xj are mapped into ti = 0, t2 = — 1> ■ ■ ■ ,tn = — ("- — !)• The 
quadrature rule for / is related to the rule for g by Qn{f) = hQn{g), with 

n 

Qn{g) = big{0) + 62ff(-l) + ■ • • + bng{l -n) = ^bk g{l - k). (14) 

fe=i 

As the integral I(q„) = /q qn{t)dt, where g„(t) = (/9„_i(t)(t + (n - 1)), is 
obviously positive, then the degree of Qn{g) (and so the degree of Qn{f)) is 
d = n — 1. Therefore, the error of Qnid) is En{g) = h ^'^"i g^"^(^) and so 

i?n(/) = /^"+'c„/(«)(e), withc„ = :?^^, eG(a,6), (15) 

n! 

assuming / to be of class C^{[a, b]). 

The Adams-Moulton rules 

Likewise, we intend to approximate /(/) = f^ f{x)dx by the rule Qn{f) = 
Yl7=i O'ifi^i)^ but now the nodes are xi = t + h, X2 = t, X3 = t — h,Xi = 
T — 2/1, • • • ,Xn = T — {n — 2)h. Performing the same change of variables as in 
the previous rule, we transform the interval [2 — n, 1] into the interval [0, 1], 
and the nodes Xi into the nodes ti = 1, ^2 = 0, . . . , t„ = 2 — n. As before 
Qnif) = hQnig). Here the polynomial qnit) is qn(t) = ^n{t){t - (2 - n)) 
and so, by Lemma [3TT] -(iii). we obtain I{qn{t)) = /q qn{t)dt / 0. Therefore 
deg(Qn(/)) = (n — 1) and from Theorem 12.11 it follows an expression for 
£"„(/) similar to (fT5D . 

Remark 3.1. In the Adams-Bashforth-Moulton rules the nodes U are in- 
tegers and so the application of the UCM leads to a linear system Ax = b, 
where ^ is a matrix of integer entries. The integrals I{ip{t)) are rational 
numbers. Therefore both the weights of Wn and the respective error coef- 
ficient Cn in ([Is]) can be computed exactly if one uses a computer algebra 
system able to represent exactly rational numbers. Of course the same is 
true for the Newton-Cotes rules discussed previously. 
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3.2 Gaussian quadrature rules 

Gaussian rules, first introduced by Gauss in [6J, rely on the properties of 
orthogonal sets of polynomials (see Gautschi [8] for a recent account on the 
subject). A Gaussian rule Gn{f) = Yl^=i'^kf{xk) is used to approximate 
the integral /(/) = /^ w{x)f{x)dx, where the so-called weight function w is 
assumed to be a positive function defined in [a, 6], continuous, and integrable 
in (a, 6). 

In the context of Gaussian quadratures the following LP' inner product 
plays a central role 

{g,h) = / w{x)g{x)h{x)dx. (16) 

J a 

The existence of an orthogonal polynomial basis {Oq,Oi, . . . , On} for P„ is 
guaranteed by the Gram-Schmit process, and it is common to call the ele- 
ments of such bases by orthogonal polynomials. The most widely used sets 
of orthogonal polynomials are known by the names of Hermite, Laguerre, 
Jacobi, Chebyshev and Legendre polynomials. Lists of orthogonal polyno- 
mials, respective weight functions and integration interval can be found, for 
instance, in Gaustchi [8], Evans [5] and Krylov [TO] , 

Gaussian quadratures are based on the choice of the nodes xi, X2, • • • , x„ 
as roots of an orthogonal polynomial 0„ of degree n. Here we do not address 
the computation of these roots, we assume they are available. 

The roots of an orthogonal polynomial are simple, real, and lie in the 
interval (a, h) (see Atkinson \i\ or Gautschi [8]). Also, as it is easy to prove, 
any orthogonal polynomial Oj of degree j is orthogonal to any polynomial 
p of degree less than j. 

The Theorem 12 . 1 1 applied to a n-point Gaussian rule enables us to recover 
well known results for the degree and error expression of this type of rules 
as shown in the next Proposition. 

Proposition 3.1. Let xi, . . . ,Xn be the roots of an orthogonal polynomial On 
of degree n. Consider L{f) = f^ w{x)f{x)dx, where w is a weight function 
and f eC^''{[a,b]). 

A Gaussian quadrature rule Gnif) = X^iLi '^«/(^«) f^''" approximating 
L{f) has degree of precision d = 2n — 1. The respective error En{f) is 

where a is the coefficient of x"^ in On, C some point in {a,b), and the norm 
\\On\\ is the L^ norm induced by the inner product (I16|) . 
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Proof. Let Bq and B be the bases, respectively for Pn-i and P2n, defined 
by ([3]) and (j3]). By Theorem 12.11 it is sufficient to prove that I{q2n) 7^ 
and I{qj) = for n < j < 2n — 1. It is obvious that I{q2n) 7^ since the 
polynomial Qn is the positive function q2n = {x — xi)'^{x — X2)^ ■ ■ ■ {x — x„)^, 
and so I{q2n) = / w{x)q2n{x)dx > 0. Hence, the degree of a n-point rule is 
at most 2n — 1. 

Let us now prove that /(g„+fc) = 0, for all < A; < n — 1. For, consider 
the following basis of P2n-i: 

B = {Oo,Oi,. . . ,On,Onfl,Onf2,- ■ ■ , On(Pn-l} , 

where {Oq, Oi, . . . , O^} is an orthogonal polynomial basis of P„. Any poly- 
nomial qn+k of degree n + k (with < A; < n— 1) can be written as a (unique) 
linear combination of the elements of B, that is 

n k 

Qn+kix) = ^7iOi(x) + ^7„+jO„(a;)(/?j(x). (18) 

For any integer < k < n — 1, the integral I{qn+k) is 
I (qn+k) = / w{x)qn+k{x)dx = {qn+k,'Po) 



n k 

i=0 j=l 

k 
= Jo{Oo,ipo) +^-/n+j{0„,ipj) =7o(Oo,(/?o), (19) 

where the last two equalities follow from the fact that an orthogonal poly- 
nomial Oj of degree j is orthogonal to any polynomial of degree less than 

3- 

Evaluating (|18p at the nodes xi, . . . ,x„, and taking into account these 
nodes are roots of both 0„ and qn+ki we obtain a homogeneous linear system 
Bx = 0, in the unknowns 70,71, .. . ,7n-i) such that the first column of B 
has all entries equal to 1. Thus, 70 = and from (fT9]l we get I (qn+k) = 0. 

It remains to show that the error expression has the form given in (J17p . 
As xi, . . . ,Xn are the roots of On, we have 0„ = a{x — xi) ■ ■ ■ {x — Xn), where 
a is the coefficient of x" in On- Then, 

rb 
|2 _ /r^ /™^ /o f^w 2 / „../^\r^2/\J 2 



||0„(a;)|| = {Onix),On{x)) = a / w{x)On{x)dx = a / w{x)q2n{x)dx. 

J a J a 

So, I{q2n) = ^\\On{x)f and (HZD follows from ©. D 
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Example: Gauss-Legendre rule 

The Legendre polynomials, usually denoted by P„, are orthogonal with re- 
spect to the weight function w{x) = 1 in the integration interval [—1,1]. 
The coefficient a and ||-Pn|P in the error expression (|17|) satisfy the follow- 
ing equalities (see, for instance Davis [3], p. 33): 

a-^^^ WPf-^ (20) 

""2"(n!)2' ll^'^ll "2n + l- ^^^> 

Thus, from (I17p . the error expression for the Gauss-Legendre rule is 

i^C,(/) = C„/<->«) "* C„ = i^;^^-JpL_, «€(-l.l). (21) 

agreeing with the well known expression for the error of a n-point Gauss- 
Legendre rule (see, for instance, Atkinson [1], p. 276). 

Computational remarks 

The computational efficiency of the algorithm referred in Remark 12.11 -2). 
either for the rules previously discussed or any other polynomial quadrature 
rule one may construct, deserve further studies. However, we have yet devel- 
oped a Mathematica (Wolfram [16]) code for the referred algorithm which, 
for a given number n of nodes, produces a list {nodes, weights, Cn}, where 
On is the coefficient in the respective error expression En{f). This code 
has been implemented for all the rules detailed in the previous section and 
the computed coefficients c„ agree with those in (fTSl) , (fT5]) and ([2T]) . When- 
ever possible the elements of the mentioned list have been computed exactly, 
such as in the cases of Newton-Cotes, Adams-Bashford and Adams-Moulton 
rules. Our program has been also tested for the Gauss-Legendre rule for a 
wide range of nodes number. For this rule, when 2 < n < 5, we obtained 
closed expressions for the respective pairs (Wn,Cn), while for 6 < n < 256 
we have computed 100-digits of precision approximations for the weights. 
The computations were carried out on a small laptop. The control of the 
error in the computation of the weights has been monitored taking into ac- 
count that the sum of the weights satisfy the equality Y27=i ^i ~ ^' ^^"^ 
knowing that all the weights in Wn are positive. Our numerical results were 
compared with those available in some published tables: for 2 < n < 512, 
Stroud and Secrest ([I5], p. 99), [BOS']; for 2 < n < 48, Krylov (fl^, p. 337), 
[205"]; for 2 < n < 128, Evans (0, p. 303), [305]; for 2 < n < 16, Kythe 
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and Schiiferkotter ([H], p. 505), [165] - the symbol kS means k significant 
digits in the corresponding table. 

We consider our high precision in situ computation of nodes, weights, 
and error coefficients for the n-point Gauss-Legendre's rule a tribute to the 
pioneering work of Lowan, Davids and Levenson [T^] and [13j . who managed 
to publish a 16-digits precision Gauss-Legendre table, for 2 < n < 16, in the 
terrible years of WWII. 

Ackno^vledgments 

MMG is supported by Instituto de Mecanica-IDMEC/IST, Centro de Pro- 
jecto Mecanico, through FCT (Portugal) /program POCTI. MESD is sup- 
ported by the FCT through the Program POCI 2010/FEDER. 

References 

[1] K. E. Atkinson, An introduction to numerical analysis. Second edition. 
John Wiley & Sons, Inc., New York, 1989. 

[2] I. S. Berezin and N. P. Zhidkov, Computing methods. Volume I, Per- 
ganrion Press, Oxford, 1965. 

[3] W. Cheney and D. Kincaid, Numerical mathematics and computing. 
Fifth edition. Brooks/Cole-Thomson, 2004. 

[4] P. J. Davis and P. Rabinowitz, Methods of numerical integration. Sec- 
ond edition. Computer Science and Applied Mathematics. Academic 
Press, Inc., Orlando, FL, 1984. 

[5] G. Evans, Practical numerical integration. John Wiley & Sons, Ltd., 
Chichester, 1993. 

[6] C. F. Gauss, Methodus nova integralium valores per approximationen 
inveniendi, C. F. Werke, Gottingen : Koniglichen Gesellschaft der Wis- 
senschaften 3, 163-196, 1886. 

[7] W. Gautschi, Numerical analysis. An introduction. Birkhduser Boston, 
Inc., Boston, MA, 1997. 

[8] W. Gautschi, Orthogonal polynomials: computation and approxima- 
tion. Numerical Mathematics and Scientific Computation. Oxford Sci- 
ence Publications. Oxford University Press, New York, 2004. 



15 



[9] P. Henrici, Discrete variable methods in ordinary differential equations. 
John Wiley & Sons, Inc., New York-London, 1962. 

[10] V. I. Krylov, Approximate Calculation of Integrals. Dover, New York, 
2005. 

[11] P. K. Kythe and M. R. Sachaferkotter, Handbook of computational 
methods for integration. Chapman & Hall/CRC, Boca Raton, FL, 2005. 

[12] A. N. Lowan, N. Davids and A. Levenson, Table of the zeros of the Leg- 
endre polynomials of order 1-16 and the weight coefficients for Gauss' 
mechanical quadrature formula. Bull. Amer. Math. Soc. 48 (1942), 739- 
743. 

[13] A. N. Lowan, N. Davids and A. Levenson, Errata to "Table of the zeros 
of the Legendre polynomials of order 1-16 and the weight coefficients 
for Gauss' mechanical quadrature formula." Bull. Amer. Math. Soc. 49 
(1943), 939. 

[14] J. F. Steffensen, Interpolation. Second Edition. Dover, Boston, 2006. 

[15] A. H. Stroud and D. Secrest, Gaussian quadrature formulas. Prentice- 
Hall, Inc., Englewood Cliffs, N.J., 1966. 

[16] S. Wolfram, The Mathematica® book. Fourth edition. Wolfram Media, 
Inc., Champaign, IL; Cambridge University Press, Cambridge, 1999. 



16 



